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Abstract 

Inhomogeneous boson systems, such as the dilute gases of integral spin atoms in low- 
temperature magnetic traps, are believed to be well described by the Gross-Pitaevskii 
equation (GPE). GPE is a nonlinear Schrodinger equation which describes the order 
parameter of such systems at the mean field level. In the present work, we describe a 
Fortran 90 computer program developed by us, which solves the GPE using a basis 
set expansion technique. In this technique, the condensate wave function (order 
parameter) is expanded in terms of the solutions of the simple-harmonic oscillator 
(SHO) characterizing the atomic trap. Additionally, the same approach is also used 
to solve the problems in which the trap is weakly anharmonic, and the anharmonic 
potential can be expressed in a polynomial in the position operators x, y, and z. The 
resulting eigenvalue problem is solved iteratively using either the self-consistent-field 
(SCF) approach, or the imaginary time steepest-descent (SD) approach. Iterations 
can be initiated using either the simple-harmonic-oscillator ground state solution, 
or the Thomas-Fermi (TF) solution. It is found that for condensates containing 
up to a few hundred atoms, both approaches lead to rapid convergence. However, 
in the strong interaction limit of condensates containing thousands of atoms, it 
is the SD approach coupled with the TF starting orbitals, which leads to quick 
convergence. Our results for harmonic traps are also compared with those published 
by other authors using different numerical approaches, and excellent agreement is 
obtained. GPE is also solved for a few anharmonic potentials, and the influence of 
anharmonicity on the condensate is discussed. Additionally, the notion of Shannon 
entropy for the condensate wave function is defined and studied as a function of 
the number of particles in the trap. It is demonstrated numerically that the entropy 
increases with the particle number in a monotonic way. 
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Program Summary 

Title of program: bose.x 
Catalogue Identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, 
N. Ireland 

Distribution format: tar.gz 

Computers : PC's/Linux, Sun Ultra 10/Solaris, HP Alpha/Tru64, IBM/AIX 
Programming language used: mostly Fortran 90 

Number of bytes in distributed program, including test data, etc.: size of the 
tar file 153600 bytes 

Number of lines in distributed program, including test data, etc. : lines in the 
tar file 4221 

Card punching code: ASCII 

Nature of physical problem: It is widely believed that the static properties of 
dilute Bose condensates, as obtained in atomic traps, can be described to a 
fairly good accuracy by the time-independent Gross-Pitaevskii equation. This 
program presents an efficient approach of solving this equation. 
Method of Solution: The solutions of the Gross-Pitaevskii equation correspond- 
ing to the condensates in atomic traps are expanded as linear combinations 
of simple-harmonic oscillator eigenfunctions. Thus, the Gross-Pitaevskii equa- 
tion which is a second-order nonlinear differential equation, is transformed 
into a matrix eigenvalue problem. Thereby, its solutions are obtained in a self- 
consistent manner, using methods of computational linear algebra. 
Unusual features of the program: None 



1 Introduction 

Ever since the discovery of Bose-Einstein condensation (BEC) in dilute atomic 
gases[l,2,3], theoretical studies of this and related phenomenon in such sys- 
tems have grown exponentially [4]. For most of the theoretical studies of BEC 
in dilute gases, the starting point is the so-called Gross-Pitaevskii equation 
(GPE)[5,6], which is nothing but a mean-field Schrodinger equation for a 
system of Bosons interacting through a two-body interaction described by 
5-function. In all but the simplest of the cases, one needs to solve the GPE 
using numerical methods. For problems involving the static properties of the 
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condensate, the numerical solutions of the time-independent GPE are of in- 
terest. And, indeed, over last several years, a significant amount of work has 
been performed towards developing novel approaches and algorithms meant for 
solving both time-dependent and independent GPE. Next we survey some of 
the recent literature in the field, restricting ourselves to the methods aimed at 
solving the time-independent GPE, which is the subject of the present paper. 
Edwards and Burnett developed a Runge-Kutta method based finite difference 
approach for solving the time-independent GPE for spherical condensates [7]. 
In another paper Edwards et al. used the basis set approach similar to the one 
presented here, to solve the GPE for anisotropic traps[8]. Dalfovo and Stringari 
developed a finite-difference based method for solving the time-independent 
GPE both for the ground state, and the vortex states, in anisotropic traps[9]. 
Esry used a finite-element approach to solve for the both the time-independent 
GPE, as well as, the Hartree-Fock equations for bosons confined in anisotropic 
traps[10]. Schneider and Feder used a discrete variable representation (DVR), 
coupled with a Gaussian quadrature integration scheme, to obtain the ground 
and the excited states of GPE in three dimensions[ll]. Adhikari used a finite- 
difference based approach to solve the two-dimensional time-independent GPE[12,13]. 
Tosi and coworkers developed finite-difference, and imaginary-time, approach 
for solving the time-independent GPE[14]. Recently, Bao and Tang developed 
a novel scheme for obtaining the ground state of the GPE, by directly minimiz- 
ing the corresponding energy functional[lo]. Additionally, utilizing harmonic 
oscillator basis set, and Gauss-Hermite quadrature integration scheme, Dion 
and Cances have proposed a scheme for solving both the time-dependent and 
-independent GPE[16]. Earlier, we had also proposed an alternative scheme for 
dealing with condensates with a large number of particles, and high number 
densities[17]. 

In this paper, we describe a Fortran program developed by us which solves 
time-independent GPE corresponding to bosons trapped in Harmonic traps. 
Instead of using the more common finite-difference approach, we have chosen 
the basis-set based approach popular in quantum chemistry [18]. The basis set 
chosen for this case is the Cartesian simple-harmonic-oscillator (SHO) basis 
set. The choice of a Cartesian basis set allows us to treat the cases ranging 
from spherical condensates to completely anisotropic condensates on an equal 
footing. Additionally, using the same approach our program allows to solve 
the time-independent GPE for anharmonic traps as well provided the anhar- 
monic term can be expressed as a polynomial in various powers of coordinates 
X, y, and z. As far as the SCF solution of the GPE is concerned, our program 
allows both the matrix-diagonalization based scheme, as well as the use of the 
imaginary-time steepest-descent method. Our program also allows the user to 
initiate the SCF process using either the SHO ground state orbital, or the 
Thomas-Fermi solution. We present the results of the calculations performed 
with our code for several interesting cases, and very good agreement is ob- 
tained with the existing results in the literature. Additionally, in the present 



3 



paper we have defined the notion of Shannon entropy in the context of GPE, 
and presented various quantitative calculations of the quantity. 



Remainder of the paper is organized as follows. In the next section we discuss 
the basic theoretical aspects of our approach. Next, in section 3, we briefly de- 
scribe the most important subroutines that comprise our program. In section 
4 we provide detailed explanation about installing and compiling our pro- 
gram. Additionally, in the same section we explain how to prepare an input 
file, and describe the contents of a typical output file. In section 5 we discuss 
the convergence properties of the program with respect to the: (a) size of the 
basis set, and (b) the method of solution. In section 6, we discuss results of 
several example runs of our program, and compare them to those published 
earlier. Additionally, in the same section, we present our results on the Shan- 
non entropy of the condensate, and on the solutions obtained in the presence 
of various anharmonic potentials. Finally, in section 7, we present our con- 
clusions. In the Appendix we present the derivation of an analytical formula 
which we have used in our program to compute the two-particle interaction 
integrals. 



2 Theory 



For the present case, the time-independent Gross-Pitaevskii equation is 

+ ^^-t W + ^^|^(r)r)^(r) = /.^(r), (1) 

where ipi^r) is condensate wave function one is solving for. VQ^^ {r) = ^m{ujlx^+ 
ujyli^ + ujIz^) + V"'"'^{x, y, z) is the confining potential for a general anisotropic 
trap {uJi^ are the trap frequencies) with the anharmonic term V^'^^{x,y,z), a 
is the s-wave scattering length characterizing the two-body interactions among 
the atoms, is the total number of bosons in the condensate, and /i is the 
chemical potential. We are assuming that the condensate wave function is nor- 
malized to unity. Before attempting numerical solutions of Eq. (1), we cast it 
in a dimensionless form by making the transformations[lo] 



f = — (2a) 



= aH^^Piv) (2b) 
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where a-c = \ is the "harmonic oscillator length" in the x direction. This 
finally leads to the dimensionless form of the GPE 



(-^v2 + v^^^{i) + K\^{mm = mi), (3) 

where u is the chemical potential in the units of fiWx, n = is a dimen- 
sionless constant determining the strength of the two-body interactions in 
the gas, and for the harmonic oscillator potential V^xt(^) ~ + lyV"^ + 
7^5^) + V'^^'^{x,y,z)^ with 7j/ = ^ and Iz = ^ being the two aspect ratios. 
We will now discuss the basis-set expansion technique, used for solving the 
GPE[8,11,16]. In this approach, one expands -(/^(f) as a linear combination of 
basis functions of three-dimensional anisotropic simple harmonic oscillator 

i=l i=l 

where (f>n^i{x), (pnyiiv), and (pn^iiz), are the harmonic oscillator basis functions 
corresponding to x, y, and z directions, respectively, Cj is the expansion coef- 
ficient, and Nbasis is the total number of basis functions used. In dimensionless 
units, e.g., 4>n^i{z) can be written as 

where H^^^^z^J^) is a Hermite polynomial of order n^j in the variable z.^J^. 
The form of the basis functions (l>n^i{x) and (bnyiiy) can be easily deduced from 
Eq. (5). Upon substituting Eq. (4) in Eq. (3), then multiplying both sides with 
another basis function ^j{x,y,z) and integrating with respect to x, y, and z 
the time-independent GPE is converted into an eigenvalue problem[18] 

FC = jlC, (6) 

where C represents the column vector containing expansion coefficients Cj's 
as its components, and the elements of the matrix F are given by 

F,j = EAj + V^r'' + 9 E Aj,k,iCkCi. (7) 

k,l=l 

Above ^ ^ ^ 

Ei - (n^i + 2) + (^yi + 2)^/y + (^^^ + g)^^' 

V^"'j'^ are the matrix elements of the anharmonic term in the confining poten- 
tial, and Jij,k,i is the boson-boson repulsion matrix. For anharmonic potentials 
which can be written as polynomials in x, y, and z, V^'^^ can be computed 
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quite easily, while Ji.j,k,i can be written as a product of three submatrices 
corresponding to the three Cartesian directions[8] 

Ji,j,k,l Jnj:ina:jnxknxiJnyinyjnyi,nyiJn^in^jn^knzl- (9) 

It can be shown that the elements of submatrices J can be written in the form 

/oo 
d^<t>nm<t>n,mn,mnm (10) 
-oo 

where (/)„;(^)'s are the harmonic oscillator basis functions of Eq. (5). Inte- 
grals involved in Eq. (10) can be computed numerically using the methods of 
Gaussian quadrature[ll,16], or analytically[8] using the formulas derived by 
Busbridge[19]. In the present work, we have used this analytical expression — 
derived in the appendix for the sake of completeness — to compute the values 
of J integrals. 

The eigenvalue problem of Eq. (6) has to be solved selfconsistently. In our pro- 
gram, for fixed values of N , this equation can be solved for fi and C using 
either the iterative diagonalization common in quantum chemistry [18], or the 
steepest-descent approach as used by Dalfovo and Stringari[9]. In both the 
approaches, the onset of self-consistency is signalled once the energy per par- 
ticle of the condensate converges to within a user-specified threshold. This 
approach is different from that of Edwards et a/. [8] where they solved Eq. (6) 
for N , using fixed values of fi. 

For the case of relatively small particle number, i.e., for a weakly interacting 
system, it does not matter what is the nature of starting guess for the con- 
densate orbital for initiating the SCF cycles. However, for the case of systems 
with large particle number, the convergence obtained is very slow (if at all), in 
case the starting condensate is taken to be the ground state of the harmonic 
trap. In such cases, the convergence is easily obtained if the starting guess for 
the condensate is taken to be of the Thomas- Fermi form, obtained by setting 
the kinetic energy term in Eq. 3 to zero 

1 / ~9 I 9 ~9 I 9 ~9\ 

= + (11) 

9 

where the Thomas-Fermi chemical potential {Itf is given by 

- ^ V a. ) ■ ^^^^ 



In our program, we can also compute the Shannon entropy associated with the 
condensate. The Shannon mixing entropy, for a general ensemble, is defined 
as[20] 

S ^ -Y^Pi^og Pi , (13) 
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where Pi is the probability for a system to be in the i-th state. Of course, in 
the present case, we do not have a thermodynamic system in which the mixing 
of various states will take place due to thermal fluctuations driven by its finite 
temperature. Thus, the question is as to how to define Shannon entropy for the 
present system, which is essentially being treated as a zero-temperature quan- 
tum mechanical system. For the purpose we adopt an information theoretic 
point-of-view, and define the probability Pi as 

Pi = \Ci\^. (14) 

where Cj are the expansion coeflacients of various SHO eigenstates in the ex- 
pression for the condensate wave function of Eq. (4). In this picture, ground 
state of the condensate is seen as a statistical mixture of the various eigen- 
states of SHO, with the mixing probability Pj. It is important to realize that 
the reason behind this mixing of states is the inter-particle interaction in the 
condensate, because, in its absence, the condensate will be in the ground state 
of the SHO (Cj = 5i^i) leading to S — 0. Thus, in a sense, entropy defined 
as per Eqs. (13) and (14) is a measure of inter-particle interactions in the 
system. Because, stronger the inter-particle interactions, the condensate will 
be a mixture of larger number of states, leading to a larger entropy. From an 
information-theoretic point-of-view larger entropy implies loss of information 
about the system, because, in such a case, the system is a mixture of a larger 
number of states. The point to be remembered, however, is that this informa- 
tion loss is being driven by the inter-particle interactions in the system while 
the corresponding information loss in a thermodynamic system is driven by 
its finite temperature, and the thermal fiuctuations caused by it. In various 
contexts, other authors have also computed and discussed the information 
entropy associated with interacting quantum systems[21,22,23]. 



3 Description of the program 

In this section we briefly describe the main program and various subroutines 
which constitute the entire module. All the subroutines, except for the di- 
agonalization subroutine taken from EISPACK[24], have been written in the 
Fortran 90 language. 

3.1 Main Program OSCL 

The main program is called OSCL, and its task is to read all the input infor- 
mation, and, among other things, fix the dimensions of various arrays. All the 
arrays needed in the program are dynamically allocated either in the main pro- 
gram, or in some of the subroutines. By utilizing the dynamic array allocation 
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facility of the Fortran 90 language, we have made the program independent of 
the size of the calculations undertaken. Because of this, the program needs to 
be compiled only once, and will run until the point the memory available on 
the computer is exhausted. Besides reading all the necessary input, program 
OSCL calls various subroutines in which different tasks associated with the 
calculation are performed. 



3.2 XMAT_0 

XMAT_0 is a small subroutine whose job is to compute the matrix elements 
of the position operator in harmonic oscillator units, with respect to the basis 
set of a one-dimensional SHO. This routine is called from the main program 
OSCL, and results are stored in a two-dimensional array called xmatrx. 



3.3 Basis Set Generation 

The calculations are performed using a basis set of a three-dimensional SHO 
consistent with the symmetry of the system. The basis set to be used is 
generated by calling one of the following three routines: (a) for a spherical 
condensate (complete isotropy) routine BASGEN3D_IS0 is used, (b) for a 
cylindrical condensate routine BASGEX3D_CYL is called, and (c) the ba- 
sis set for a completely anisotropic condensate is generated using the routine 
BASGEN3D_ANIS0. In all the cases the basis functions are arranged in the 
ascending order of their harmonic oscillator energies and, if needed, a heap 
sort routine called HPSORT is used to achieve that. All these subroutines have 
the option of imposing parity symmetry on the basis set if the potential has 
that symmetry. This leads to a substantial reduction in the size of the basis 
set in most cases. 



3.4 HAM0_3D: 

This subroutine is called from the main program OSCL and its purpose is to 
generate the matrix elements of the noninteracting (one-particle) part of the 
condensate Hamiltonian. If the condensate is confined in a perfectly harmonic 
trap, one-particle part of the Hamiltonian is trivial. However, for the case 
where the trap potential is anharmonic, the potential matrix elements are 
generated from the position operator matrix elements xmatrx{i,j) mentioned 
above. 
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3.5 BEC_DRV 



Subroutine BEC DRV is called from the main program OSCL, and as its 
name suggests, it is the driver routine for performing calculations of the con- 
densate using the time-independent GPE. Apart from allocating a few arrays, 
the main task of this routine is to call either: (a) routine BOSE_SCF meant for 
solving for the condensate wave function using the iterative-diagonalization- 
based SCF approach, or (b) routine BOSE_STEEP used for solving for the 
condensate wave function using the steepest-descent approach of Dalfovo and 
Stringari[9]. 

3.6 BOSESCF 

This subroutine solves the time-independent GPE in a self-consistent manner 
using an iterative diagonalization approach. Its main tasks are as follows: 

(1) Allocate various arrays needed for the SCF calculations 

(2) Setup the starting orbitals. For this the options are: (i) diagonalize the 
one-particle part of the Hamiltonian, (ii) use the Thomas- Fermi orbitals, 
or (iii) use the orbitals obtained in a previous run. 

(3) Perform the SCF calculations. For the purpose, the two-particle inte- 
grals Ji,j,k,i (cf. Eq. (9) are calculated on the fly during each iteration 
using the formulas derived in the appendix. In other words the storage of 
these matrix elements is completely avoided, thereby saving substantial 
amount of memory and disk space. This approach is akin to the "direct 
SCF" approach utilized in quantum chemistry. Permutation symmetries 
of indices i, j, /c, and / are utilized to reduce the number of integrals 
evaluated. Moreover, the evaluation of an integral is undertaken only if 
it is found to be nonzero as per the symmetry selection rules. The in- 
tegrals in question are evaluated in a subroutine called JMNPQ_CAL. 
The F matrix constructed in each iteration is diagonalized through a 
Householder diagonalization routine called HOUSEH, which is from the 
EISPACK package of routines[24], and is written in Fortran 77. 

(4) The chemical potential and the condensate wave function obtained after 
every iteration are written in various data files so that the progress of the 
calculation can be monitored. 

3.7 BOSESTEEP 

Alternatively, the condensate wave function and the chemical potential can be 
obtained using the subroutine BOSE_STEEP which, instead of the iterative 
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diagonalization approach, utilizes the steepest-descent approach to achieve 
convergence, starting from a given starting orbital. In this approach, as out- 
lined by Dalfovo and Stringari[9], the starting orbitals are evolved towards 
the true orbitals in small imaginary time steps by repeated application of the 
Hamiltonian, i.e., the F operator. Here, the main computational step is the 
multiplication of a column vector by a matrix, which in the present version of 
the program is achieved by a call to the Fortran 90 intrinsic function MAT- 
MUL. However, one could certainly try to improve upon this by developing a 
subroutine which can utilize the symmetric nature of the Fock matrix. Apart 
from this, rest of the actions performed in this subroutine are identical to 
those of BOSE_SCF. 

3.8 THOMAS_FERMI 

This subroutine is invoked either from the subroutine BOSE_SCF or from 
BOSE_STEEP in case the SCF calculations are to be initiated by assuming 
Thomas-Fermi form of the starting orbitals. Upon invocation, this subrou- 
tine directly constructs the operator F corresponding to the Thomas-Fermi 
orbitals. In this case, the r-space integration is performed using a trapezoidal- 
rule-based scheme on a three-dimensional Cartesian grid. 

3.9 ENTROPY 

This subroutine is called if the entropy of the condensate needs to be computed 
with respect to the Harmonic oscillator basis functions, as per Eqs. 13 and 14. 
It is a very small subroutine with a straightforward implementation. 

3.10 CONDPLOT 

This subroutine computes the numerical values of the condensate wave func- 
tion for a user-specified set of points in space. The numerical values of the 
Hermite polynomials needed for the purpose are computed using the subrou- 
tine HERMITE, described below. 

3.11 HERMITE 

This subroutine computes the values of the Hermite polynomial Hn{x) for a 
set of user specified values of x, and order n. Fast computation of polynomials 
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is achieved by using initializations Hq{x) — 1, Hi{x) — 2x, and the recursion 
relation Hn+i{x) — 2xHn{x) — 2nHn-i{x) . 



4 Installation 

All the files needed to install and run the program are kept in the gzipped, 
tarred archive bose.tar.gz. It consists of: (a) All the Fortran files containing 
the main program (file oscl .f 90), and various subroutines, called by the main 
program, (b) four versions of Makefiles which can be used for compiling the 
code on various Linux/Unix systems, and (c) several sample input and output 
files in a subdirectory called Examples. The program was developed on Pen- 
tium 4 based machine running Redhat Fedora core 1 operating system using 
noncommercial version of the Intel Fortran compiler version 8.1. However, it 
has also been verified that it runs on Sun Solaris Sparc based systems, Compaq 
alpha (now HP alpha) based systems running True Unix, and IBM Power PC 
systems running AIX. For these systems, the Fortran 90 compilers supplied 
with those operating systems were used. In order to install and compile the 
program, following steps need to be followed: 

(1) Uncompress the program files in a directory of user's choice using the 
command gunzip bose .tar .gz followed by tar -xvf bose.tar. 

(2) Verify that the four makefiles Makef ile_linux, Makef ile_solaris, Makef ile_alpha, 
and Makef ile_aix are present. Copy the suitable version of the make file 

to the file Makefile. For example, if the system is a Sun Solaris Sparc 
system, issue the command cp Makef ile_solaris Makefile. 

(3) Now issue the command make which will initiate the compilation. If ev- 
erything is successful, upon completion bin directory of your account will 
have the program execution file bose.x. If your account does not have 
a directory named bin, you will have to either create this directory, or 
modify the Makefile to ensure that the file bose.x is created in the 
directory of your choice. 

(4) If the bin directory is in your path, try running the program using one of 
the sample input files located in the subdirectory Examples. For exam- 
ple, by issuing the command bose.x < bec_iso.dat > x.out one can 
run the program for an isotropic trap and the output will be written 
in a file called x.out. This should be compared with the supplied file 
bec_iso.out to make sure that results obtained agree with those of the 
example run. 

Additionally, a file called README is also provided which lists and briefly ex- 
plains all the files included in the package. Although we have not investigated 
the installation of the program on operating systems other than Linux/Unix, 
we do not anticipate any problems with such operating systems. 
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4-1 Input Files 

In order to keep the input process as free of errors as possible, we have adopted 
the philosophy that before each important input cards, there will be a com- 
pulsory comment line. It is irrelevant as to what is written in the comment 
lines, but, by writing something meaningful, one can keep the input process 
transparent. The input quantities following the comment line have to be in free 
format, with the restriction that the ASCII input cards should be in uppercase 
letters. Because of the use of comment lines, the input files are more or less 
self explanatory. In the sample input files, we have started all the comment 
lines with the character #. A sample input file corresponding to a cylindrical 
trap potential is listed below 

#Type of oscillator 

CYLINDRICAL 

# NXMAX, NYMAX, NZMAX 
10, 8 

# NO. OF TERMS IN THE ANHARMONIC POTENTIAL 


# OMEGAX, OMEGAY, OMEGAZ, JILA parameters 
1.0, 2.8284271 

# Type of SCF equation 
GP 

# No . of particles 
1000 

# Scattering Length, JILA Rb87 
4.33d-3 

# SCF convergence threshold, Maximum # of allowed iterations. 
l.d-8, 1000 

# Whether Parity is a good quantum number or not 
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PARITY 



# Method for calculations 
SCF 

# Starting orbitals 
SHO 

# Whether orbital Mixing will be done 
FOCKMIX 

0.4 

# Whether orbital plots needed 
PLOT 

-5.0,5.0,0.05 



1,0,0 

# Entropy Calculation 
ENTROPY 



Next we describe the input cards one by one. 

(1) First card is an ASCII card describing the type of trap potential. Options 
are: ISOTROPIC, CYLINDRICAL, or ANISOTROPIC. 

(2) Second card specifies the maximum quantum numbers of the basis func- 
tions to be included for various directions. For an isotropic trap one entry 
is needed {n.j. = Uy = n^), for cylindrical trap, two entries are needed 
{rix = Uy and n^), while for an anisotropic oscillator three entries are 
needed {hx, Uy, Uz). These numbers eventually determine the total num- 
ber of basis functions N^asis used to expand the condensate as per Eq. 
4. 

(3) Third card deals with the anharmonic terms in the trap potential. Any 
potential of the form J2^i'' Cja;™^'!/™*'^;™^ can be added to the Harmonic 
trap potential. The first entry here is Nanh after which Nanh entries con- 
sisting of {Ci,m^ ,m^} are given. In the example input, no anhar- 
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monicity was considered, thus Nanh has been set to zero. 

(4) Fourth card deals with the trap frequencies with the convention that 
ujx = ^, and rest of the frequencies measured in the units of uUx- The 
example input corresponds to the trap frequencies of the JILA experiment 
with oja^^ojy^l, uj^^ 2.284271. 

(5) Fifth card specifies which mean-field equation is to be solved. Options 
are GP for the Gross-Pitaevskii equation, and HF for the Hartree-Fock 
equation. 

(6) Sixth card reads the total number of bosons in the trap. 

(7) Seventh card is the value of the s-wave scattering length in the Harmonic 
oscillator units. 

(8) Eighth card inputs the convergence threshold, followed by the maximum 
number of iterations allowed to achieve convergence. 

(9) Ninth card specifies whether parity should be treated as a good quantum 
number or not. Options are PARITY and NOPARITY. If the trap poten- 
tial is invariant under the parity operation, use of this card leads to a 
tremendous reduction in the size of the basis set needed to solve for the 
condensate wave function. 

(10) Tenth card specifies as to which method is to be used for solving the 
mean-field equations. Options are SCF corresponding to the iterative di- 
agonalization method, and STEEPEST-DESCENT corresponding to the use 
steepest-descent method of Dalfovo and Stringari[9]. In case one opts for 
the steepest-descent approach, the size of the time step to be used in the 
calculations also needs to be specified. 

(11) Eleventh card specifies as to what sort of starting guess for the con- 
densate should be used to start the solution process. Valid options are 
SHO corresponding to the simple-harmonic oscillator ground state and 
THOMAS-FERMI corresponding to the Thomas-Fermi form of the conden- 
sate. 

(12) It has been found that in several difficult cases, convergence can be 
achieved if one utilizes the techniques of Fock matrix mixing or conden- 
sate orbital mixing[9,10]. Valid options are (a) FOCKMIX for Fock matrix 
mixing (b) ORBMIX for condensate mixing, (c) Any other ASCII entry 
such as NOMIX for neither of these options. In case options (a) or (b) are 
chosen, one needs to specify the parameter xmix quantifying the mixing 
according to the formula 

R^'^ = xmix R^^ + (1 - xmix) R^'-^\ 

where is the quantity under consideration in the i-th iteration. Thus, 
if Fock matrix mixing has been opted, xmix specifies the fraction of the 
new Fock matrix in the total Fock matrix in the i-th iteration. 

(13) This is the penultimate card which decides whether the user wants the 
numerical values of condensate wave function along user specified set of 
data points, such that the condensate could be plotted as a function of 
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spatial coordinates. Keyword PLOT means that the answer is in affirma- 
tive while any other option such as NOPLOT will disable the numerical 
computation of the condensate. If the keyword PLOT has been supplied 
as in the example input, further data rmin, rmax, dr is specified next 
which determines the starting position, ending position, and the step size 
for generating the points on which the condensate is to be computed. 
After these values, we need to specify variable ndir which is the number 
of directions along which the condensate needs to be computed. Finally, 
ndir Cartesian directions have to be specified. The example input file 
instructs the program to compute the value of the condensate along the 
X axis, for —5.0 < a; < 5.0, in the steps of 0.05. 
(14) Final card specifies whether one wants to compute the mixing entropy. 
Valid options are ENTROPY and any other entry such as NOENTROPY. If 
one opts for entropy calculation, one can do so for a whole range of 
eigenfunctions specified by their lower bound and the upper bound. Entry 
1 , 1 specifies that entropy of only the ground state needs to be computed. 



4-2 Output file 

Apart from the usual information related to various system parameters, the 
most important information that an output file contains is the approach (or 
lack thereof) to convergence of the calculations as far as the chemical potential 
is concerned. Besides that, any other computed quantity such as the entropy 
is also listed in the output file. The important portions of the output file, cor- 
responding to the input file discussed in the previous section, are reproduced 
below. The complete sample output file is called bec_cyl_jila. out, and is 
included in the tar archive. 

SCF iterations begin 



Starting chemical potential= 2.4142135 



Iteration # 


Chem. Pot. 


Energy/particle 


Energy-Converg 


1 


3.876643 


5.3193762 


5.3193762 


2 


4.271632 


4.0689872 


-1.2503890 


3 


4.517259 


3.8580062 


-0.2109810 


4 


4.619476 


3.8439807 


-0.0140255 


5 


4.689213 


3.8418453 


-0.0021354 



15 



6 



4.720534 



3.8413295 



-0.0005159 



4.743576 



3.8411720 



-0.0001575 



4.753914 



3.8411212 



-0.0000508 



4.761791 



3.8411040 



-0.0000172 



10 



4.765292 



3.8410982 



-0.0000059 



11 



4.768022 



3.8410961 



-0.0000020 



12 



4.769223 



3.8410954 



-0.0000007 



13 



4.770177 



3.8410951 



-0.0000003 



14 



4.770592 



3.8410950 



-0.0000001 



15 



4.770927 



3.8410950 



0.0000000 



16 



4.771071 



3.8410950 



0.0000000 



17 



4.771189 



3.8410950 



0.0000000 



Convergence achieved on the BEC ground state 

Eigenstate # Information Entropy 

1 0.7477159 

The contents of the output file listed above are self-explanatory. It basically 
shows that after seventeen iterations, the total energy per particle of the con- 
densate converges to the value 3.841095, leading to the chemical potential 
value of 4.771. Additionally, the entropy of the condensate is computed to be 
0.744771. 

In addition to the above mentioned main output file, there is another output 
file created in the ASCII format which contains the condensate orbital ob- 
tained at the end of each iteration. This file is written in the logical unit 9, 
and is named orbitals .dat. When a new run is started, the program always 
looks for this file and tries to use the condensate solution present there to start 
the iterations. In other words, condensate solution present in orbitals.dat 
is used to restart an old aborted run. If some incompatibility is found between 
the condensate solution, and the present run, the solution in the orbitals . dat 
is ignored and a new run is initiated. Thus, if one wants to start a completely 
new run, any old orbitals.dat file must first be deleted. 
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5 Convergence Issues 

In this section we compare the convergence of our results with respect to 
the size of the basis set used. We also compare the convergence properties of 
different iterative approaches aimed at obtaining the condensate solutions. 



5.1 Convergence with respect to the basis set 

Before treating the results obtained as the true results, one must be sure as to 
convergence properties with respect to the size of the basis set. This aspect of 
the calculations is explored in the present section by means of two examples 
corresponding to condensates in spherical and cylindrical traps, respectively, 
with fifteen hundred {N — 1500) bosons each. 

First we discuss the case of the condensate in an isotropic trap, results for 
which are presented in table 1. The value of the scattering length used in the 
calculations is listed in the caption of the table. For this case, only one value 
specifying the largest quantum number nmax for the basis functions, needs to 
be specified. It is obvious from the table that the results which have converged 
to three decimal places both in the chemical potential and the entropy require 
nmax = 8, leading to the total number of basis functions Nhasis = 35. This 
means that the size of the Fock matrix diagonalized during the iterative di- 
agonalization is 35 x 35, which is computationally very inexpensive. It is also 
obvious from the table that in order to get four decimal place convergence, we 
only need to use nmax = 10 corresponding to a 56 x 56 Fock matrix, whose 
diagonalization can also be carried out quite fast. 

Similar results for the condensate in a cylindrical trap corresponding to the 
JILA parameters[2o] are presented in table 2. Because of the anisotropy of 
the trap, the convergence is to be judged with respect to two parameters 
nxmax deciding the highest quantum number of the basis functions for x— 
and y— directions, and nzmax the corresponding number for the 2;-direction. 
We will first try to understand the convergence properties using a few heuristic 
arguments. In the JILA experiment [25] the trap frequency in the 2;-direction 
ujz was more than twice the value of the trap frequencies in the x- and y- 
directions, ujx and ujy. Therefore, due to inter-particle repulsion, the condensate 
will be much more delocalized along the x/y-directions, as compared to the z- 
direction. This means that, in order to achieve convergence, one would expect 
to use higher-energy basis functions in the x/y-directions, as compared to the 
;i;-direction. In other words, at convergence nxmax > nzmax. And when we 
examine table 2, we find that this is indeed the case. We notice that the three- 
decimal place convergence in both the chemical potential, and the entropy. 
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is obtained for nxmax = 8 and nzmax — 6, although the data presented in 
the table covers a much larger range of parameters. Thus, we conclude that 
reasonably accurate values of various physical quantities can be obtained with 
basis sets of modest sizes, both for the isotropic as well as for the cylindrical 
condensates. 



5.2 Comparison of different numerical approaches 

As mentioned earlier, our program can solve the GP equation using two numer- 
ical approaches: (a) iterative diagonalization (ID) of the Fock matrix, and (b) 
steepest-descent (SD) method of Dalfovo et al.{9]. In the previous sections, all 
the presented results were obtained by the ID method. In the present section, 
we would like to present results based upon the steepest-descent approach, and 
compare them to those obtained using the iterative diagonalization method. 
We present our results for the cylindrical trap corresponding to the JILA 
parameters[25], with an increasing number of particles in table 3. Obviously, 
the numerical solution of the GP equation becomes increasingly difficult as the 
number of particles in the condensate grows, because of the increased contri- 
bution of the inter-particle repulsion. Therefore, it is very important to know 
as to how various numerical approaches perform as N is gradually increased. 

As the results presented in the table suggest that for smaller values of N, 
neither of the two approaches have any problems achieving convergence, and 
the results obtained were found to be in agreement with each other to three 
decimal places for the basis set used. We found that in most of the cases, 
the ID approach worked only when Fock-matrix mixing used. Although, we 
managed to achieve convergence for the cases depicted in table 3 with the ID 
method; however, as the number of bosons in the condensate grows further, 
the convergence becomes slow and difficult to achieve by this method, a fact 
also emphasized by Schneider and Feder[ll]. On the other hand, the SD ap- 
proach did not have any convergence problems for the cases we investigated. 
As depicted in table 3, the SD approach led to convergence both with the 
SHO starting orbital, as well as the Thomas-Fermi starting orbital. However, 
for larger values of A'", the use of Thomas-Fermi solution as the starting guess 
for the condensate, will lead to much faster convergence with this approach. 
Thus, we conclude that: (a) For smaller values of A^, both the ID as well as the 
SD methods will lead to convergence, and (b) for really large values of A^, the 
convergence is guaranteed only with the SD method. In the SD method the 
main computational operation is the multiplication of a vector by a matrix, 
which will be significantly faster as compared to the matrix diagonalization 
procedure needed by the ID method for calculations involving large basis sets. 
Thus, for calculations involving large basis sets, SD method will be faster as 
compared to the ID method. Therefore, all things considered, we believe that 
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the SD method is the more robust of the two possible approaches. 



6 Example Runs 

In this section we report the results of our calculations for a variety of trap 
parameters, and compare our results with those published by other authors. 
We also study the behavior of the Shannon entropy of the condensate with 
respect to the number of particles it contains. Additionally, we also present 
our results for the cases of anharmonic traps. 

6.1 Comparison with other works 

In this section we present the results of several calculations performed on both 
isotropic and the cylindrical traps, and compare them to the results obtained 
by other authors. Several authors have performed such calculations, however, 
for the sake of brevity, we restrict our comparisons mainly to the works of Bao 
and Tang[15] for the spherical condensate, and to the results of Dalfovo and 
Stringari[9] for the cylindrical condensate. 

Recently, using finite-element based approach, Bao and Tang[lo] performed 
calculations for condensates on a variety of harmonic traps, and presented re- 
sults as a function of the interaction parameter k — ^ttoTV table 4 our results 

. . ."^ 

for the chemical potentials of condensates in isotropic traps corresponding to 
increasing values of k are compared with those reported by Bao and Tang[lo]. 
The agreement between the results obtained by two approaches is exact to the 
decimal places reported by Bao and Tang[15]. Note that the aforesaid agree- 
ment was obtained for rather modest basis set sizes, and calculations were 
completed on a personal computer in a matter of minutes. 

Next we discuss the results obtained for a cylindrical trap corresponding to 
JILA parameters[25] for an increasing number of bosons. Our results are pre- 
sented in table 5, where they are also compared to the results of Dalfovo and 
Stringari obtained using a finite-difference based approach [9]. The agreement 
between our results and those of Dalfovo and Stringari is virtually exact for 
all their reported calculations[9]. Again, the noteworthy point is that this level 
of agreement was obtained with the use of modest sized basis sets, and the 
computer time running into a few minutes. 

Thus, the excellent agreement between our results, and those obtained by 
other authors using different approaches, gives us confidence about the essen- 
tial correctness of our methodology. Now the question arises, will this numer- 
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ical method work for values of the interaction parameter k which are much 
larger than the ones considered here. The encouraging aspect of the approach 
is that for none of the larger values of n which we considered did we experi- 
ence a numerical breakdown of the approach. It is just that for larger values 
of K, the total number of basis functions needed to achieve convergence on the 
chemical potential will be larger as compared to the smaller This, of 

course, will also lead to an increase in the CPU time needed to perform con- 
verged calculations. For example, for the case of the spherical trap considered 
in table 4, when we doubled k to the value 6274, we needed to use basis func- 
tions corresponding to nmax — 20 with Nj,asis — 286 to achieve two-decimal 
places convergence in the chemical potential. For k, = 9411, to achieve similar 
convergence, these numbers increased to nmax = 24 and Nbasis = 455. Finally, 
when K, was increased to 15685, the corresponding numbers were nmax = 28 
and Nbasis = 680, with the CPU time running into several hours. For the 
cylindrical trap (cf. table 5), for N — 20000 bosons {k — 1088.2) the conver- 
gence was achieved in a matter of minutes with nxmax — 14, nzmax = 8 
and Nbasis = 180. When the number of bosons in the trap was doubled such 
that K = 2176.4, similar level of convergence on the chemical potential was 
obtained with nxmax = 18, nzmax = 8 and Nbasis = 275. Even with much 
larger values of k (> 30000) both for the spherical, and the cylindrical traps, 
we did not encounter any convergence difficulties when the calculations were 
performed with the modest sized basis sets mentioned earlier. But it was quite 
obvious that, to obtain highly accurate values of chemical potentials for such 
cases, one will have to use basis sets running into thousands which will make 
the calculations quite time consuming. 

At this point, we would also like to compare our approach to that of Schneider 
and Feder[ll], who used a DVR based technique to obtain accurate solutions 
of the time-independent GPE. In the DVR approach the basis functions are 
the so-called "coordinate eigenfunctions", which, in turn, are assumed to be 
linear-combinations of other functions such as the SHO eigenstates, or the 
Lagrange interpolating functions[ll]. Thus in the DVR approach of Schnei- 
der and Feder[ll], the SHO eigenstates are used as intermediate basis func- 
tions, and not as primary basis functions as is done in our approach. Us- 
ing this approach, coupled with the "direct-inversion in the iterative space" 
(DIIS) method, Feder and Schneider managed to obtain accurate solutions 
for anisotropic condensates for quite large values of the interaction parameter 
K,{11]. However, the price to be paid for this accuracy was the use of a very 
large basis set consisting of several thousands of basis functions[ll] even for 
rather small values of k. 

Finally, we present the plots of the condensates in a spherical trap, for increas- 
ing values of A^, in Fig. 1. As expected, the calculations predict a depletion 
of central condensate density, and corresponding delocalization of the con- 
densate, with increasing N. The results presented are in excellent agreement 
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with similar results presented by various other authors[9,15]. Moreover, if we 
compare the value of the condensate at the center of the trap (|V'(0, 0, 0)|) for 
the isotropic trap with the published results of Bao and Tang[15], we again 
obtain excellent agreement for all values of A^. 



6.2 Anharmonic Potentials 

Recently, several studies have appeared in the literature studying the influence 
of trap anharmonicities on the condensates, in light of rotating condensates, 
and the resultant vortex structure[26]. However, we approach the influence 
of trap anharmonicity from a different perspective, namely that of quantum 
chaos. Therefore, the anharmonicities considered here are in the absence of 
any rotation, and the aim is to study their influence on the ground and the 
excited states of the condensate. We assume the unperturbed harmonic trap to 
be the cylindrical one corresponding to the JILA parameters[25], and consider 
two types of anharmonic perturbations in the x ~ y plane: (a) the Henon- 
Heiles potential with V"^(a;, j/) = a{x^y — ^y^), and (b) the Fourleg poten- 
tial V°'"''^(x,y) — ax'^y'^, where a is the anharmonicity parameter. Note that 
the Henon-Heiles potential reduces the circular symmetry of the cylindrical 
trap in the x — y plane to the triangular one (symmetry group C^^), and 
the fourleg potential reduces the symmetry to that of a square (group C41,). 
In case of Henon-Heiles potential the inversion symmetry of the cylindrical 
trap is also destroyed, while for the fourleg potential, it is still preserved. The 
Henon-Heiles potential introduces deconflnement in the trap, the Fourleg po- 
tential, on the other hand, strengthens the confinement of the original trap. 
Both these potentials are known to lead to chaotic behavior for higher energy 
states, both at the classical and quantum-mechanical levels of theories[27,28]. 
In a separate work communicated elsewhere, we have examined the excited 
states of condensates under the influence of these potentials, in order to an- 
alyze the signatures of chaotic behavior. In the present work, however, we 
intend only to demonstrate the capabilities of our program as far as the an- 
harmonicity is concerned, and restrict ourselves only to the ground states of 
the condensates in presence of these potentials. Results of our calculations on 
the chemical potentials of condensates in a cylindrical trap corresponding to 
the JILA experiment [25], and = 1000, are presented in table 6 as a function 
of anharmonicity a. Corresponding plots of the condensate along the y axis are 
presented in Fig. 2. As far as the influence of anharmonicity on the chemical 
potential is concerned, from table 6 we conclude that, for a given value of N, 
for increasing a, it increases for the Henon-Heiles potential, and decreases for 
the fourleg potential. Similarly, upon examining the Fig. 2, we conclude that 
for the Henon-Heiles potential, the central condensate density gets depleted 
with increasing a, while the behavior in case of the fourleg potential is just the 
opposite. Additionally, the fact that the inversion symmetry is broken in case 



21 



of Henon-Heiles potential, is obvious from the asymmetry of corresponding 
condensate plots. 



6.3 Entropy Calculations 



Here we discuss the Shannon entropy of condensates in isotropic and cylin- 
drical traps, as a function of the dimensionless strength parameter k = 
Since the scattering length in most of the traps is fixed, for such cases the 
change in k can be construed as due to changes in N. In Fig. 3 we present the 
plots of Shannon entropy versus k plots for the condensates trapped both in 
isotropic, as well as cylindrical traps. Although a detailed analysis of the Shan- 
non entropy of condensates is being presented elsewhere, we make a couple of 
important observations: (i) For both types of traps the entropy increases as a 
function of k. Initially, the rate of increase is quite high, but for larger values 
of K, it settles down to a much lower value, (ii) For a given nonzero value of k, 
the entropy of a condensate in a cylindrical trap is always larger than that of 
a condensate in a spherical trap. In other words, the trap anisotropy appears 
to increase the Shannon entropy of the system. 



7 Conclusions 



In this paper we have reported a Fortran 90 implementation of a harmonic 
oscillator basis set based approach towards obtaining the numerical solutions 
of time independent GPE. We have presented applications of our program to 
a variety of situations including anharmonic potentials, and in calculations of 
the Shannon entropy of the condensate. We also compared the results obtained 
from our program to those obtained by other authors, and found near-perfect 
agreement. Therefore, we encourage the users to apply our program to a variety 
of situations, and contact us in case they encounter errors. We have extensive 
plans for further development of our program. Some of the possible directions 
are: (a) extension of our approach to time-dependent GPE, allowing one to 
deal with condensate dynamics, (b) taking condensate rotation into account, 
allowing one to study the vortex phenomena, and (c) dealing with condensates 
with nonzero spins, i.e., the so-called spinor condensates[29]. Work along these 
lines is presently in progress in our group, and, upon completion, will be 
reported in future publications. 
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A Appendix 



Here our aim is to compute the two-particle integral of Eq. (10) defined as 

/oo 
-oo 

where, in terms of the dimensionless coordinates ^, the single particle wave 
function 4>ni{0 given by 



MO 



(A.2) 



Substituting Eq.A.2 in Eq. (10), we get 



JfiinjUkni 



'■ninjUkTii -I 



(A.3) 



4,n,n,n.= / e-^^'i/^, (O^/n, (O^n, (0^n,(Orfe- (A.4) 



where 



Since Hermite polynomials have a definite parity, the integral ImnjUkni will be 
nonvanishing only if the sum + Uj + n/,. + nj is an even number. Now we will 
use a standard result for the product of two Hermite Polynomials, 



min{m,ra} 
k=0 







ln\ 









where 



m 
k 



etc. are the binomial coefficients. Upon substituting Eq. A. 5 in 



Eq. A.4, we obtain 

min{m,?i} min{i3,r} 



'■m,n,q,r 



m \ I n \ I q 



E E 2'=^'^!^! ■ MM 
fe=o 1=0 \ k I \k J \l 



f \ roo 



/oo 
-oo 

(A.6) 

In order to perform the integral above, we recall the result derived by Busbridge[19] 



e-^« Hm{i)Hn{i)di = (-1) — 2^r( ), (A.7) 
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for m + n = even. Substituting this we get 



m+rt— p — g m+n+p+q— 1 



-1)^-2 



m,n,q,ri 



where 





fn] 




(A 




(;) 


[k) 


[k) 




V) 







min{m,n} min{g,r} 
fe=0 1=0 



Upon substituting Eqs. (A. 8) and (A. 9) in Eq. (A. 3), we get 

1 



(A.8) 



,171 + n + q + r — 2k — 21 1 



J, 



m,n,q,r 



m+n—q—r 
= (-1) 2 



Tr^2{m\n\q\r\) 



m,n,q,r- 



2 

(A.9) 
(A.IO) 



On using the expression 



r(^) 



2n! 
22"n! 



(A.ll) 



and setting ■m + n + q + r = 2t, we have 

m + n + 5 + r-2A;-2Z 1 {2t - 2k - 21)1^^ 



+ 7;) 



2' 22«-2fc-2'(i- A;-/)!' 



(A.12) 



Upon substituting Eq. (A.12) and the values of binomial coefficients in Eq. 
(A.9), we have 



n 



m,n,q,r 2m+n+q+r 



E 



m\n\q\rlx 
(^-iy-k2^{k+i)^2t-2k-2l)\ 



(A.13) 



y{m- ky.{n - ky.{q - iy.{r - iy{k\){u){t -k-iy: 

which, upon substitution in Eq. (A.IO), leads to the final expression 



J, 



(— 1) ^2'' lm\n\q\r\ 



m,n,q,r 2m+n+q+r 



27r 



-'m,n,q,r-i 



(A.14) 



where 



(-iy-k2^k+2i(^2t-2k-2iy 



min{m,n} nim{q,r} 

S S (^A;)!(n-fc)!(g-/)!(r-/)!(fc!)(/!)(i-A;-Or 

(A.15) 



Expressions of Eqs. (A.14) and (A.15) have been used in the function JINTT, 
which is called via subroutine JMNPQ_CAL, to compute these two-body in- 
tegrals. We would like to emphasize that the series of Eq. (A.15) has terms 
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with alternating signs, and, therefore, is potentially unstable for large values of 
m, n, q, and r. Thus, it is crucial to use high arithmetic precision while sum- 
ming the series. With the usual double-precision arithmetic (REAL*8 vari- 
ables), we found that the the series was unstable for values of m, n, q, r 
larger than 16. To circumvent these problems, we used quadruple precision 
(REAL* 16 variables) in function JINTT to sum the series. Once the summa- 
tion is performed, the results are converted into the double-precision format. 
We believe that this approach has made the two-particle integral calculation 
process very robust, and accurate. 
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Table 1 

Convergence on the chemical potential and the mixing entropy of the condensate 
in an isotropic trap with scattering length a = 2.4964249 x 10~^ax and number of 
bosons N = 1500, with respect to the basis set size, nmax is the maximum value of 
the quantum number of the SHO basis function in a given direction, Nhasis is the 
total number of basis functions corresponding to a given value of nmax, and N^er 
represents the total number of SCF iterations needed to achieve convergence on 
the condensate energy per particle. In all the calculations iterative diagonalization 
method, along with Fock matrix mixing with xmax = 0.6, was used. The SCF 
convergence threshold was 1.0 x 10~^. 



nmax 






Chemical Potential 


Entropy 


2 


4 


19 


2.939116 


0.5083669 


4 


10 


11 


2.915046 


0.5387074 


6 


20 


14 


2.911181 


0.5396975 


8 


35 


14 


2.911278 


0.539088 


10 


56 


15 


2.911375 


0.5392423 


12 


84 


15 


2.911346 


0.5392770 


14 


120 


15 


2.911337 


0.5392822 


16 


165 


13 


2.911337 


0.5392797 
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Table 2 

Convergence on the chemical potential and the mixing entropy of the condensate in 
a cylindrical trap with trap parameters corresponding to the JILA experiment [25], 
and the number of bosons N = 1500, with respect to the basis set size, nxmax is 
the maximum value of the quantum number of the SHO basis ftmction in x- and 
direction, nzmax is the same number corresponding to the z-direction. Rest of 
the quantities have the same meaning as explained in the caption of table 1. In 
all the calculations iterative diagonalization, along with Fock matrix mixing with 
xmax = 0.3, was used. The SCF convergence threshold was 1.0 x 10~ ^. 



nxmax 


nzmax 


-^basis 


-^iter 


Chemical Potential 


Entropy 


4 





6 


39 


5.786323 


0.9387124 


4 


2 


12 


19 


5.421930 


0.9383227 


4 


4 


18 


27 


5.423981 


0.9380013 


6 





10 


31 


5.737602 


0.9770727 


6 


2 


20 


19 


5.405221 


0.9549194 


6 


4 


30 


20 


5.405440 


0.9545441 


6 


6 


40 


20 


5.404970 


0.9545920 


8 





15 


21 


5.737074 


0.9752491 


8 


2 


30 


20 


5.404605 


0.9546530 


8 


4 


45 


21 


5.404691 


0.9542773 


8 


6 


60 


21 


5.404218 


0.9543291 


8 


8 


75 


21 


5.404147 


0.9543239 


10 





21 


21 


5.736935 


0.9755365 


10 


2 


42 


20 


5.404741 


0.9544814 


10 


4 


63 


20 


5.404540 


0.9541044 


10 


6 


84 


20 


5.404066 


0.9541561 


10 


8 


105 


18 


5.403994 


0.9541508 


10 


10 


126 


18 


5.403991 


0.9541477 


12 





28 


21 


5.736763 


0.9756539 


12 


2 


56 


20 


5.404636 


0.9545789 


12 


4 


84 


20 


5.404436 


0.9542004 


12 


6 


112 


20 


5.403962 


0.9542517 


12 


8 


140 


20 


5.403891 


0.9542463 


12 


10 


168 


20 


5.403888 


0.9542432 


12 


12 


196 


18 


5.403888 


0.9542428 
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Table 3 

Comparison of the chemical potentials (jj) obtained using the iterative diagonaliza- 
tion (ID) technique, and the steepest-descent (SD) technique[9], for the condensates 
in a cylindrical trap with trap parameters corresponding to the JILA experiment [25], 
and a given number of bosons (N). Quantities Nn^r, nxmax, and nzmax have the 
same meaning as in the caption of table 2, and n is expressed in the units of hux- 
In the ID based calculations for N < 2000, SHO ground state solution was used to 
start the iterations, while for larger values of N the iterations were started using 
the Thomas-Fermi approximation. In all the cases corresponding to the ID method, 
Fock matrix mixing was used, with 0.05 < xmxix < 0.3. In the SD based calcula- 
tions, the iterations were started using the Tlionias-Fermi approximation, with the 
size of the time step being 0.02 units. 



N 


nxmax 


nzmax 




NiteriSB) 


m(id) 


m(sd) 


500 


8 


6 


16 


62 


3.938611 


3.938865 


1000 


8 


6 


14 


88 


4.770707 


4.772055 


1500 


8 


6 


21 


95 


5.404218 


5.405491 


2000 


12 


8 


23 


104 


5.931870 


5.932878 


10000 


14 


8 


86 


99 


10.505267 


10.505124 


15000 


14 


8 


105 


129 


12.239700 


12.239465 


20000 


14 


8 


104 


109 


13.665923 


13.665686 



Table 4 

Comparison of the chemical potentials (in the imits of hw) obtained from our pro- 
gram, and those reported by Bao and Tang[15], for an isotropic trap, with increasing 
values of interaction parameter k. The negative value of k implies attractive inter- 
particle interactions. For the value of scattering length stated in table 1, k = 3137.1 
corresponds to A'^ = 1 X 10^ bosons. Symbols nmax and Ajjer have the same mean- 
ing as in the previous tables. For the last two calculations, SD method with a time 
step of 0.02 units, and Thomas- Fermi initial guess were employed. Our chemical 
potentials have been truncated to as many decimal places as reported by Bao and 
Tang[15l. 



K 


nmax 




/x(This work) 


/x(Ref.[15]) 


-3.1371 


14 


120 


1.2652 


1.2652 


3.1371 


14 


120 


1.6774 


1.6774 


12.5484 


14 


120 


2.0650 


2.0650 


31.371 


14 


120 


2.5861 


2.5861 


125.484 


14 


120 


4.0141 


4.0141 


627.42 


16 


165 


7.2485 


7.2484 


3137.1 


16 


165 


13.553 


13.553 
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Table 5 

Comparison of the chemical potentials (in the units of hwx) obtained from our pro- 
gram, and those reported by Dalfovo and Stringari[9], for a cylindrical trap corre- 
sponding to the JILA parameters[25], with increasing number A'^ of bosons. Symbols 
nxmax, nzmax, and Niter have the same meaning as in the previous tables. Cal- 
culations for N > 10000 were performed by the SD method using Thomas-Fermi 
starting orbitals, a time-step of 0.02 units, and a convergence threshold of 1.0 x 10^'''. 
We have truncated our chemical potentials to as many decimal places as reported 
by Dalfovo and Stringari[9]. 



N 


nxmax 


nzmax 


Nhasis 


/i(This work) 


;u(Ref.[9]) 


100 


8 


6 


60 


2.88 


2.88 


200 


8 


6 


60 


3.21 


3.21 


500 


8 


6 


60 


3.94 


3.94 


1000 


8 


6 


60 


4.77 


4.77 


2000 


8 


6 


60 


5.93 


5.93 


5000 


10 


8 


105 


8.14 


8.14 


10000 


10 


8 


105 


10.5 


10.5 


15000 


14 


8 


180 


12.2 


12.2 


20000 


14 


8 


180 


13.7 


13.7 



Table 6 

Influence of trap anharmonicities on the chemical potential. The table below presents 
results for the Henon-Heiles, and the fourleg potentials, for cylindrical trap corre- 
sponding to JILA parameters[25], with A'^ = 1000. 



a 


/i(Henon-Heiles) 


/Lt (Fourleg) 


0.00 


4.7712 


4.7712 


0.03 


4.7662 


4.8261 


0.06 


4.7497 


4.8752 


0.09 


4.7207 


4.9202 


0.12 


4.6764 


4.9619 


0.15 


4.6131 


5.0009 
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Figure 1. Plots of condensates along the x axis for an isotropic trap, with an in- 
creasing number of N of bosons. The trap parameters used were the same as in the 
data of tables 1 and 4. Lines correspond to N = 100, 500, 1000, 5000, and 10000, 
and are in the descending order of the central condensate density, and distances are 
in the units of a-r- 



Fig 1 : Tiwari and Shukia 



0.4 - 



0.3 



0.2 - 



0.1 - 







-5 




-4 -3 -2 




X 



12 3 4 



31 



Figure 2. Influence of various types of anharmonicities on condensates in cylin- 
drical traps with A'^ = 1000, and scattering length corresponding to the JILA 
parameters [25]. The plots correspond to: (a) the Henon-Heiles potential, and (b) 
the Fourleg potential. In each graph, solid, dotted, and dashed lines represent values 
of anharmonicity parameter (see text) a = 0.0, 0.05, and 0.15, respectively. The y 
coordinate is measured in the units of ax- 

Fig 2: Tiwari and Shukia 
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Figure 3. Plots of Shannon entropy of condensates in an isotropic trap (solid line) and 
cylindrical trap (dashed line), as a function of the dimensionless strength parameter 

^ _ iTTaN 
— — I • 

Fig 3: Tiwari and Shukia 
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